load packages

library(tidyverse)
library(knitr)

define palettes

pal_control = c("#2E86B4", "#F2B827")
pal_outcome = wesanderson::wes_palette("Zissou1", 3, "continuous")
pal_condition = wesanderson::wes_palette("Zissou1", 6, "continuous")
pal_condition8 = c("#3B9AB2", "#357382", "#E86900", "#7FB8BA", "#ED4100", "#CC1400", "#EECC20", "#E1B002")
pal_condition8 = c("#CC1400", "#ED4100", "#E86900", "#357382", "#3B9AB2", "#7FB8BA", "#E1B002", "#EECC20")

source data

  • Note that % fat variable was incorrect in 2012_FDES –> replaced it with the data from the EMA dataframe
source("load_data.R")
saveRDS(ind_diffs, "ind_diffs.RDS")
saveRDS(ema, "ema.RDS")

ind_diffs = ind_diffs %>%
  mutate(bmi_z = as.numeric(scale(bmi)),
         fat_z = as.numeric(scale(fat)))

standardize brain data

  • Windsorizing observations +/- 3 SD from the mean to 3 –> otherwise we lose these subs in the SCA analyses
  • Z-score within ROI and session for ROIs
  • Z-score within map, test (association or uniformity), mask (masked or unmasked) and session for PEVs
# roi distributions
betas %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, roi, process, meanPE) %>%
  unique() %>%
  filter(!is.na(roi)) %>%
  ggplot(aes(meanPE, fill = roi)) +
  geom_density(color = NA) +
  facet_wrap(~roi, ncol = 4) +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none")

# dot product distribution
dots %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = test)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(.85, .2))

dots %>%
  filter(session == "all") %>%
  select(subjectID, con, condition, control, process, test, dotProduct, mask) %>%
  unique() %>%
  ggplot(aes(dotProduct, fill = mask, alpha = .5)) +
  geom_density(color = NA) +
  facet_wrap(~process, ncol = 3, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = c(.85, .2))

# standardize and winsorize
betas_std = betas %>%
  group_by(roi, session) %>%
  mutate(meanPE_std = scale(meanPE, center = TRUE, scale = TRUE),
         outlier = ifelse(meanPE_std > 3, "yes",
                   ifelse(meanPE_std < -3, "yes", "no")),
         meanPE_std = ifelse(meanPE_std > 3, 3,
                      ifelse(meanPE_std < -3, -3, meanPE_std))) %>%
  ungroup()

dots_std = dots %>%
  group_by(map, test, mask, session) %>%
  mutate(dotProduct_std = scale(dotProduct, center = TRUE, scale = TRUE),
         outlier = ifelse(dotProduct_std > 3, "yes",
                   ifelse(dotProduct_std < -3, "yes", "no")),
         dotProduct_std = ifelse(dotProduct_std > 3, 3,
                          ifelse(dotProduct_std < -3, -3, dotProduct_std))) %>%
  ungroup()

summarize outliers

# summarize
betas_std %>%
  filter(session == "all") %>%
  group_by(outlier) %>%
  summarize(n = n()) %>%
  spread(outlier, n) %>%
  mutate(percent = round((yes / (yes + no)) * 100, 2))
dots_std %>%
  filter(session == "all") %>%
  group_by(outlier, test, mask) %>%
  summarize(n = n()) %>%
  spread(outlier, n) %>%
  mutate(percent = round((yes / (yes + no)) * 100, 2))
# remove outlier column
betas_std = betas_std %>%
  select(-outlier)

dots_std = dots_std %>%
  select(-outlier)

merge data

dataset = full_join(betas_std, dots_std, by = c("subjectID", "con", "process", "condition", "control", "session")) %>%
  filter(!condition == "social" & mask == "unmasked") %>%
  mutate(subjectID = as.character(subjectID)) %>%
  left_join(., ind_diffs, by = "subjectID") %>%
  ungroup() %>%
  mutate(subjectID = as.factor(subjectID),
         condition = as.factor(condition),
         control = as.factor(control),
         roi = as.factor(roi),
         process = as.factor(process),
         test = as.factor(test))

saveRDS(dataset, "full_dataset.RDS")
saveRDS(betas_std, "betas_dataset.RDS")
saveRDS(dots_std, "dots_dataset.RDS")

tidy data for plotting

  • Select cons generated across all sessions
  • Calculate mean across ROIs within each psychological process
  • Calculate balance score ROI = cognitive control - reward
ema_tidy = ema %>%
  select(subjectID, enact_prop) %>%
  unique() %>%
  filter(!is.na(subjectID)) %>%
  mutate(enact_prop_z = as.numeric(scale(enact_prop)))

dots_plot = dots_std %>%
  filter(!condition == "social" & mask == "unmasked" & session == "all") %>%
  select(-c(map, con, mask, session, dotProduct)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender)) %>%
  left_join(., ema_tidy, by = "subjectID")

betas_plot = betas_std %>%
  filter(!condition == "social" & session == "all") %>%
  group_by(subjectID, process, condition, control) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE)) %>%
  select(-c(con, xyz, roi, meanPE, sdPE, meanPE_std)) %>%
  unique() %>%
  spread(process, meanProcessPEstd) %>%
  gather(process, meanProcessPEstd, cognitive_control, reward, value) %>%
  select(-c(session)) %>%
  left_join(., ind_diffs) %>%
  select(-c(sample, DBIC_ID, age, gender)) %>%
  ungroup() %>%
  left_join(., ema_tidy, by = "subjectID")

combined_plot = dots_plot %>%
  rename("value_std" = dotProduct_std) %>%
  mutate(test = ifelse(process == "craving_regulation", "neural signature", test)) %>%
  select(subjectID, process, test, condition, control, value_std, fat, bmi, enact_prop) %>%
  bind_rows(betas_plot %>%
              rename("value_std" = meanProcessPEstd) %>%
              mutate(test = "ROI") %>%
              select(-contains("_z")))

tidy data for modeling

model_data = dataset %>%
  filter(test == "association" & session == "all" & control == "nature" & condition == "food") %>%
  unique() %>%
  group_by(process, subjectID) %>%
  mutate(meanProcessPEstd = mean(meanPE_std, na.rm = TRUE),
         processPE = paste0(process, "_ROI"),
         processPEV = paste0(process, "_PEV")) %>%
  ungroup() %>%
  select(-c(xyz, meanPE, meanPE_std, sdPE, dotProduct, map, process)) %>%
  spread(processPEV, dotProduct_std) %>%
  spread(processPE, meanProcessPEstd) %>%
  group_by(subjectID) %>%
  fill(contains("PEV"), .direction = "down") %>%
  fill(contains("PEV"), .direction = "up") %>%
  fill(contains("ROI"), .direction = "down") %>%
  fill(contains("ROI"), .direction = "up") %>%
  select(-c(roi, craving_ROI, craving_regulation_ROI, test, mask, age, gender)) %>%
  unique()

model_data_bmi = model_data %>%
  select(-contains("fat"), -restraint) %>%
  na.omit()

model_data_fat = model_data %>%
  select(-contains("bmi"), -restraint) %>%
  na.omit()

model_data_ema = model_data %>%
  select(-contains("fat"), -contains("bmi"), -restraint) %>%
  right_join(., ema_tidy, by = "subjectID") %>%
  na.omit()

visualize

DVs and dietary restraint

# source raincloud plot
source("https://gist.githubusercontent.com/benmarwick/2a1bb0133ff568cbe28d/raw/fb53bd97121f7f9ce947837ef1a4c65a73bffb3f/geom_flat_violin.R")

# plot
ind_diffs %>%
  left_join(., ema_tidy, by = "subjectID") %>%
  select(subjectID, bmi, fat, enact_prop, restraint) %>%
  rename("BMI" = bmi,
         "body fat percentage" = fat,
         "enactment" = enact_prop,
         "dietary restraint" = restraint) %>%
  gather(var, val, -subjectID) %>%
  mutate(var = factor(var, levels = c("BMI", "body fat percentage", "enactment", "dietary restraint"))) %>%
  ggplot(aes("", val, fill = var, color = var)) +
    geom_flat_violin(position = position_nudge(x = .1, y = 0), color = FALSE) +
    geom_boxplot(width = .1, outlier.shape = NA, color = "black", fill = "white") +
    geom_point(position = position_jitter(width = .05), size = .5, alpha = .75) +
    facet_wrap(~var, scales = "free", nrow = 1) +
    scale_fill_manual(values = c(pal_outcome, pal_condition[3])) +
    scale_color_manual(values = c(pal_outcome, pal_condition[3])) +
    scale_x_discrete(expand = c(.1, 0)) +
    coord_flip() +
    labs(y = "\nvalue", x = "") +
    theme_minimal(base_size = 14) +
    theme(legend.position = "none")

combined

combined_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature")),
         process = gsub("_", " ", process),
         test = factor(test, levels = c("neural signature", "ROI", "association", "uniformity"))) %>%
  ggplot(aes(control, value_std, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_wrap(~ test + process, nrow = 3) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

combined_plot %>%
  filter(!condition == "nature") %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature")),
         process = gsub("_", " ", process),
         process = ifelse(!test == "ROI", sprintf("%s PEV", process), sprintf("%s ROI", process)),
         process = factor(process, levels = c("reward ROI", "reward PEV", "craving PEV",
                                              "cognitive control ROI", "cognitive control PEV", "craving regulation PEV",
                                              "value ROI", "value PEV")),
         test = ifelse(!test %in% c("association", "uniformity"), "none", test),
         test = factor(test, levels = c("none", "association", "uniformity"))) %>%
  ggplot(aes(control, value_std, fill = process)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~condition) +
  scale_fill_manual(name = "", values = pal_condition8) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

ROIs

means

Mean parameter estimates as a function of process, stimulus type, and control condition

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(condition, meanProcessPEstd, fill = control)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(~process) +
  scale_fill_manual(name = "", values = pal_control) +
  labs(x = "\nstimulus", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(control, meanProcessPEstd, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(~process) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

correlations

Correlations between variables and BMI, % fat, or EMA enactment proportion

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(bmi, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nBMI", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(fat, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

betas_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(enact_prop, meanProcessPEstd, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control~process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\ndaily enactment proportion", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

PEVs

means

Mean expression as a function of process, stimulus type, control condition, and test

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(condition, dotProduct_std, fill = control)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~process) +
  scale_fill_manual(name = "", values = pal_control) +
  labs(x = "\nstimulus", y = "mean expression\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "nature"))) %>%
  ggplot(aes(control, dotProduct_std, fill = condition)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "bar", position = position_dodge(.9)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "errorbar", position = position_dodge(.9), width = 0) +
  facet_grid(test~process) +
  scale_fill_manual(name = "", values = pal_condition) +
  labs(x = "\ncontrol condition", y = "mean expression\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

correlations

Correlations between variables and BMI, % fat, or EMA enactment proportion

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(bmi, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nBMI", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(fat, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

dots_plot %>%
  mutate(condition = factor(condition, levels = c("food", "dessert", "meal", "snack", "social", "nature"))) %>%
  ggplot(aes(enact_prop, dotProduct_std, color = condition)) +
  geom_point(alpha = .2) +
  geom_smooth(method = "lm", alpha = .25) +
  facet_grid(control + test ~ process) +
  scale_color_manual(name = "", values = pal_condition) +
  labs(x = "\nbody composition (% fat)", y = "mean standardized parameter estimate\n") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top",
        legend.box = "horizontal")

correlations among ROIs and PEVs

PEVs v ROIs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID, -contains("bmi"), -contains("fat"), -contains("restraint"), -contains("enact_prop")) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(!grepl("dot", Var1) & !grepl("roi", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

ROIs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(grepl("roi", Var1) & grepl("roi", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

PEVs

dots_plot %>%
  unite(var, process, test, condition, control, sep = "_") %>%
  mutate(var = sprintf("dot_%s", var)) %>%
  spread(var, dotProduct_std) %>%
  left_join(., betas_plot) %>%
  unite(var, process, condition, control, sep = "_") %>%
  mutate(var = sprintf("roi_%s", var)) %>%
  spread(var, meanProcessPEstd) %>%
  select(-subjectID) %>%
  cor(., use = "pairwise.complete.obs") %>%
  reshape2::melt() %>%
  filter(grepl("dot", Var1) & grepl("dot", Var2)) %>%
  ggplot(aes(Var1, Var2, fill = value)) +
    geom_tile(color = "white") +
    scale_fill_gradientn(name = "correlation\n", colors = c(pal_condition[1], "white", pal_condition[6]), limits = c(-1, 1), breaks = c(-1, 0, 1)) + 
    geom_text(aes(label = round(value, 1)), size = 2) +
    labs(x = "", y = "") + 
    theme(legend.text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank(),
          axis.text = element_text(color = "black"),
          axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
          axis.line = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.border = element_blank(),
          panel.background = element_blank())

descriptives

means and sds

model_data %>%
  left_join(., model_data_ema) %>%
  gather(Variable, val, contains("PEV"), contains("ROI"), bmi, fat, enact_prop, restraint) %>%
  filter(!Variable == "balance_ROI") %>%
  group_by(Variable) %>%
  summarize(M = mean(val, na.rm = TRUE),
            SD = sd(val, na.rm = TRUE)) %>%
  filter(!is.na(Variable)) %>%
  mutate(Variable = gsub("_", " ", Variable),
         Variable = ifelse(Variable == "bmi", "BMI",
                    ifelse(Variable == "fat", "Body fat percentage",
                    ifelse(Variable == "enact prop", "Enactment",
                    ifelse(Variable == "restraint", "Dietary restraint", Hmisc::capitalize(Variable))))),
         Variable = factor(Variable, levels = c("BMI", "Body fat percentage", "Enactment",
                                                "Reward ROI", "Reward PEV", "Craving PEV",
                                                "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
                                                "Value ROI", "Value PEV",
                                                "Dietary restraint"))) %>%
  arrange(Variable) %>%
  kable(format = "pandoc", digits = 2)
Variable M SD
BMI 23.19 3.36
Body fat percentage 25.43 8.00
Enactment 0.47 0.22
Reward ROI 0.06 0.50
Reward PEV 0.07 0.77
Craving PEV -0.13 0.76
Cognitive control ROI -0.13 0.50
Cognitive control PEV -0.11 0.78
Craving regulation PEV 0.40 0.67
Value ROI 0.40 0.68
Value PEV 0.28 0.77
Dietary restraint 16.90 8.10

correlations

cor_mat = model_data %>%
  left_join(., model_data_ema) %>%
  select(subjectID, bmi, fat, enact_prop, contains("ROI"), contains("PEV"), restraint) %>%
  gather(var, val, contains("PEV"), contains("ROI"), bmi, fat, enact_prop, restraint) %>%
  mutate(var = gsub("_", " ", var),
         var = ifelse(var == "bmi", "BMI",
               ifelse(var == "fat", "Body fat percentage",
               ifelse(var == "enact prop", "Enactment",
               ifelse(var == "restraint", "Dietary restraint", Hmisc::capitalize(var))))),
         var = factor(var, levels = c("BMI", "Body fat percentage", "Enactment",
                                      "Reward ROI", "Reward PEV", "Craving PEV",
                                      "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
                                      "Value ROI", "Value PEV",
                                      "Dietary restraint"))) %>%
  arrange(var) %>%
  spread(var, val) %>%
  ungroup() %>%
  select(-subjectID) %>%
  as.matrix()

cors = Hmisc::rcorr(cor_mat)

p = cors$P %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  gather(var, p_value, -rowname) %>%
  mutate(code = ifelse(between(p_value, .05, .1), "†",
                ifelse(between(p_value, .01, .05), "*", 
                ifelse(between(p_value, .001, .01), "**",
                ifelse(p_value < .001, "***", "")))))

cors$r %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  gather(var, val, -rowname) %>%
  left_join(., p) %>%
  mutate(val = gsub("0.(.*)", ".\\1", sprintf("%.2f%s", val, code)),
         val = gsub("1..NA", "--", val),
         rowname = factor(rowname, levels = c("BMI", "Body fat percentage", "Enactment",
                                      "Reward ROI", "Reward PEV", "Craving PEV",
                                      "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
                                      "Value ROI", "Value PEV",
                                      "Dietary restraint"))) %>%
  select(-p_value, -code) %>%
  spread(var, val) %>%
  rename(" " = rowname) %>%
  select(" ", "BMI", "Body fat percentage", "Enactment",
         "Reward ROI", "Reward PEV", "Craving PEV",
         "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
         "Value ROI", "Value PEV",
         "Dietary restraint") %>%
  kable(format = "pandoc")
BMI Body fat percentage Enactment Reward ROI Reward PEV Craving PEV Cognitive control ROI Cognitive control PEV Craving regulation PEV Value ROI Value PEV Dietary restraint
BMI .66*** -.13 -.12† -.02 -.02 -.02 -.06 .10 .06 .00 .24***
Body fat percentage .66*** -.20† -.07 -.07 -.07 -.03 -.06 .02 .01 -.04 .36***
Enactment -.13 -.20† .21* .31** .21* .01 .05 .03 .23* .29** .14
Reward ROI -.12† -.07 .21* .62*** .47*** .31*** .41*** .08 .30*** .60*** -.02
Reward PEV -.02 -.07 .31** .62*** .81*** .47*** .59*** .19** .58*** .91*** .01
Craving PEV -.02 -.07 .21* .47*** .81*** .45*** .54*** .24*** .51*** .76*** .10
Cognitive control ROI -.02 -.03 .01 .31*** .47*** .45*** .91*** .21*** .06 .30*** -.15*
Cognitive control PEV -.06 -.06 .05 .41*** .59*** .54*** .91*** .23*** .16** .40*** -.13*
Craving regulation PEV .10 .02 .03 .08 .19** .24*** .21*** .23*** .28*** .15* .02
Value ROI .06 .01 .23* .30*** .58*** .51*** .06 .16** .28*** .65*** .14*
Value PEV .00 -.04 .29** .60*** .91*** .76*** .30*** .40*** .15* .65*** .07
Dietary restraint .24*** .36*** .14 -.02 .01 .10 -.15* -.13* .02 .14* .07

t-tests

model_data %>%
  select(-contains("z"), -contains("age"), -contains("restraint"), -bmi, -fat) %>%
  ungroup() %>%
  select_if(is.numeric) %>%
  map_df(~ broom::tidy(t.test(.x, mu = 0, na.action = "na.omit")), .id = "x") %>%
  rename("Variable" = x,
         "t" = statistic,
         "df" = parameter,
         "p" = p.value) %>%
  mutate(`Mdiff [95% CI]` = sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         Variable = gsub("_", " ", Variable),
         Variable = Hmisc::capitalize(Variable),
         Variable = factor(Variable, levels = c("Reward ROI", "Reward PEV", "Craving PEV",
                                                "Cognitive control ROI", "Cognitive control PEV", "Craving regulation PEV",
                                                "Value ROI", "Value PEV"))) %>%
  select(Variable, `Mdiff [95% CI]`, t, df, p) %>%
  arrange(Variable) %>%
  kable(format = "pandoc", digits = 2)
Variable Mdiff [95% CI] t df p
Reward ROI 0.06 [-0.00, 0.12] 1.90 261 .059
Reward PEV 0.07 [-0.02, 0.17] 1.51 261 .132
Craving PEV -0.13 [-0.22, -0.04] -2.74 261 .007
Cognitive control ROI -0.13 [-0.19, -0.07] -4.14 261 < .001
Cognitive control PEV -0.11 [-0.20, -0.01] -2.20 261 .029
Craving regulation PEV 0.40 [0.31, 0.48] 9.52 261 < .001
Value ROI 0.40 [0.31, 0.48] 9.38 261 < .001
Value PEV 0.28 [0.18, 0.37] 5.84 261 < .001

run cross-validated models

specify model variables

options(na.action = "na.fail")
data_ctrl = caret::trainControl(method = "repeatedcv", number = 5, repeats = 5)
set.seed(123)

BMI

ROI models

bmi_roi_full = lm(bmi_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_bmi)
(bmi_roi_mods = MuMIn::dredge(bmi_roi_full))
MuMIn::get.models(bmi_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_bmi_roi = caret::train(bmi_z ~ reward_ROI + value_ROI,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

bmi_pev_full = lm(bmi_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_bmi)
(bmi_pev_mods = MuMIn::dredge(bmi_pev_full))
MuMIn::get.models(bmi_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_bmi_pev = caret::train(bmi_z ~ craving_regulation_PEV,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

bmi_null = lm(bmi_z ~ 1, data = model_data_bmi)
summary(bmi_null)
## 
## Call:
## lm(formula = bmi_z ~ 1, data = model_data_bmi)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2100 -0.6542 -0.0558  0.4827  4.6115 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01306    0.06326   0.206    0.837
## 
## Residual standard error: 1.006 on 252 degrees of freedom
bmi_combined = caret::train(bmi_z ~ reward_ROI + value_ROI + craving_regulation_PEV,
                     data = model_data_bmi,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(bmi_aic = AIC(bmi_null, best_bmi_roi$finalModel, best_bmi_pev$finalModel, bmi_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*bmi_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
  arrange(AIC))

summarize the overall best fitting

# ROI
summary(best_bmi_roi$finalModel)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1355 -0.6606 -0.0899  0.4662  4.5052 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.03113    0.07236  -0.430   0.6674  
## reward_ROI  -0.30523    0.12990  -2.350   0.0196 *
## value_ROI    0.15583    0.09654   1.614   0.1077  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9976 on 250 degrees of freedom
## Multiple R-squared:  0.02502,    Adjusted R-squared:  0.01722 
## F-statistic: 3.207 on 2 and 250 DF,  p-value: 0.04214
best_bmi_roi
## Linear Regression 
## 
## 253 samples
##   2 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 202, 203, 202, 201, 204, 203, ... 
## Resampling results:
## 
##   RMSE       Rsquared    MAE      
##   0.9991722  0.03203531  0.7489974
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
sd(best_bmi_roi$resample$Rsquared)
## [1] 0.04890756
sd(best_bmi_roi$resample$RMSE)
## [1] 0.1078594
# combined
summary(bmi_combined$finalModel)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1994 -0.6654 -0.0886  0.5098  4.4311 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            -0.07322    0.07768  -0.942   0.3469  
## reward_ROI             -0.30708    0.12961  -2.369   0.0186 *
## value_ROI               0.11695    0.09989   1.171   0.2428  
## craving_regulation_PEV  0.14370    0.09792   1.467   0.1435  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9953 on 249 degrees of freedom
## Multiple R-squared:  0.03338,    Adjusted R-squared:  0.02173 
## F-statistic: 2.866 on 3 and 249 DF,  p-value: 0.03725
bmi_combined
## Linear Regression 
## 
## 253 samples
##   3 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 203, 201, 203, 202, 203, 203, ... 
## Resampling results:
## 
##   RMSE      Rsquared    MAE      
##   1.004276  0.04192202  0.7616355
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
sd(bmi_combined$resample$Rsquared)
## [1] 0.045541
sd(bmi_combined$resample$RMSE)
## [1] 0.1267501

% fat

ROI models

fat_roi_full = lm(fat_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_fat)
(fat_roi_mods = MuMIn::dredge(fat_roi_full))
MuMIn::get.models(fat_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_fat_roi = caret::train(fat_z ~ reward_ROI,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

fat_pev_full = lm(fat_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_fat)
(fat_pev_mods = MuMIn::dredge(fat_pev_full))
MuMIn::get.models(fat_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_fat_pev = caret::train(fat_z ~ craving_PEV,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

fat_null = lm(fat_z ~ 1, data = model_data_fat)

fat_combined = caret::train(fat_z ~ reward_ROI + craving_PEV,
                     data = model_data_fat,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(fat_aic = AIC(fat_null, best_fat_roi$finalModel, best_fat_pev$finalModel, fat_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*fat_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
    arrange(AIC))

summarize the overall best fitting

summary(fat_null)
## 
## Call:
## lm(formula = fat_z ~ 1, data = model_data_fat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9531 -0.5945  0.1396  0.6457  2.4201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01679    0.06321   0.266    0.791
## 
## Residual standard error: 0.9994 on 249 degrees of freedom

EMA enactment proportion

ROI models

enact_prop_roi_full = lm(enact_prop_z ~ value_ROI + reward_ROI + cognitive_control_ROI, data = model_data_ema)
(enact_prop_roi_mods = MuMIn::dredge(enact_prop_roi_full))
MuMIn::get.models(enact_prop_roi_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_enact_prop_roi = caret::train(enact_prop_z ~ value_ROI,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

PEV models

enact_prop_pev_full = lm(enact_prop_z ~ value_PEV + reward_PEV + cognitive_control_PEV + craving_PEV + craving_regulation_PEV, data = model_data_ema)
(enact_prop_pev_mods = MuMIn::dredge(enact_prop_pev_full))
MuMIn::get.models(enact_prop_pev_mods, subset = TRUE) %>%
  tibble() %>%
  rename("model" = ".") %>%
  mutate(tidied = purrr::map(model, broom::tidy),
         model_num = row_number()) %>%
  select(model_num, tidied) %>%
  unnest()
best_enact_prop_pev = caret::train(enact_prop_z ~ reward_PEV,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

combined

enact_prop_null = lm(enact_prop_z ~ 1, data = model_data_ema)

enact_prop_combined = caret::train(enact_prop_z ~ value_ROI + reward_PEV,
                     data = model_data_ema,
                     trControl = data_ctrl,
                     method = "lm",
                     na.action = na.omit)

(enact_aic = AIC(enact_prop_null, best_enact_prop_roi$finalModel, best_enact_prop_pev$finalModel, enact_prop_combined$finalModel) %>%
    rownames_to_column() %>%
    extract(rowname, "model", ".*prop_(.*)\\$.*") %>%
    mutate(model = toupper(model),
           model = ifelse(is.na(model), "Null", model)) %>%
  arrange(AIC))

summarize the overall best fitting

summary(best_enact_prop_pev$finalModel)
## 
## Call:
## lm(formula = .outcome ~ ., data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.14970 -0.57410 -0.01054  0.53196  2.37344 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.06291    0.09858  -0.638  0.52488   
## reward_PEV   0.36683    0.11493   3.192  0.00192 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9573 on 95 degrees of freedom
## Multiple R-squared:  0.09686,    Adjusted R-squared:  0.08735 
## F-statistic: 10.19 on 1 and 95 DF,  p-value: 0.001917
best_enact_prop_pev
## Linear Regression 
## 
## 97 samples
##  1 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 78, 77, 78, 77, 78, 78, ... 
## Resampling results:
## 
##   RMSE       Rsquared   MAE      
##   0.9565301  0.1467801  0.7619614
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
sd(best_enact_prop_pev$resample$Rsquared)
## [1] 0.1304739
sd(best_enact_prop_pev$resample$RMSE)
## [1] 0.1571684

tables

bmi

bmi_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_bmi_roi$resample$Rsquared), 
                            mean(best_bmi_pev$resample$Rsquared), 
                            mean(bmi_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_bmi_roi$resample$Rsquared), 
                               sd(best_bmi_pev$resample$Rsquared), 
                               sd(bmi_combined$resample$Rsquared)),
                     RMSE = c(mean(best_bmi_roi$resample$RMSE), 
                              mean(best_bmi_pev$resample$RMSE), 
                              mean(bmi_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_bmi_roi$resample$RMSE), 
                                 sd(best_bmi_pev$resample$RMSE), 
                                 sd(bmi_combined$resample$RMSE)))

bmi_table = broom::tidy(best_bmi_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_bmi_roi$finalModel)$fstatistic[1], 
                             summary(best_bmi_roi$finalModel)$fstatistic[2], 
                             summary(best_bmi_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_bmi_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_bmi_pev$finalModel)$fstatistic[1], 
                             summary(best_bmi_pev$finalModel)$fstatistic[2], 
                             summary(best_bmi_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(bmi_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(bmi_combined$finalModel)$fstatistic[1], 
                             summary(bmi_combined$finalModel)$fstatistic[2], 
                             summary(bmi_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(bmi_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., bmi_aic) %>%
  left_join(., bmi_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "BMI") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

bmi_table %>%
  kable(format = "pandoc", digits = 2)
outcome model AIC r2 (SD) RMSE (SD) term b [95% CI] SE t p
BMI Combined* 721.56 0.04 (0.05) 1.00 (0.13) Intercept -0.07 [-0.23, 0.08] 0.08 -0.94 .347
BMI Combined* 721.56 0.04 (0.05) 1.00 (0.13) Reward ROI -0.31 [-0.56, -0.05] 0.13 -2.37 .019
BMI Combined* 721.56 0.04 (0.05) 1.00 (0.13) Value ROI 0.12 [-0.08, 0.31] 0.10 1.17 .243
BMI Combined* 721.56 0.04 (0.05) 1.00 (0.13) Craving regulation PEV 0.14 [-0.05, 0.34] 0.10 1.47 .144
BMI ROI* 721.74 0.03 (0.05) 1.00 (0.11) Intercept -0.03 [-0.17, 0.11] 0.07 -0.43 .667
BMI ROI* 721.74 0.03 (0.05) 1.00 (0.11) Reward ROI -0.31 [-0.56, -0.05] 0.13 -2.35 .020
BMI ROI* 721.74 0.03 (0.05) 1.00 (0.11) Value ROI 0.16 [-0.03, 0.35] 0.10 1.61 .108
BMI PEV 723.45 0.02 (0.03) 0.99 (0.14) Intercept -0.05 [-0.19, 0.10] 0.07 -0.67 .507
BMI PEV 723.45 0.02 (0.03) 0.99 (0.14) Craving regulation PEV 0.16 [-0.03, 0.34] 0.09 1.64 .102
BMI Null 724.15 Intercept 0.01 [-0.11, 0.14] 0.06 0.21 .837

% fat

fat_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_fat_roi$resample$Rsquared), 
                            mean(best_fat_pev$resample$Rsquared), 
                            mean(fat_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_fat_roi$resample$Rsquared), 
                               sd(best_fat_pev$resample$Rsquared), 
                               sd(fat_combined$resample$Rsquared)),
                     RMSE = c(mean(best_fat_roi$resample$RMSE), 
                              mean(best_fat_pev$resample$RMSE), 
                              mean(fat_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_fat_roi$resample$RMSE), 
                                 sd(best_fat_pev$resample$RMSE), 
                                 sd(fat_combined$resample$RMSE)))

fat_table = broom::tidy(best_fat_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_fat_roi$finalModel)$fstatistic[1], 
                             summary(best_fat_roi$finalModel)$fstatistic[2], 
                             summary(best_fat_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_fat_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_fat_pev$finalModel)$fstatistic[1], 
                             summary(best_fat_pev$finalModel)$fstatistic[2], 
                             summary(best_fat_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(fat_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(fat_combined$finalModel)$fstatistic[1], 
                             summary(fat_combined$finalModel)$fstatistic[2], 
                             summary(fat_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(fat_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., fat_aic) %>%
  left_join(., fat_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "Body fat") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

fat_table %>%
  kable(format = "pandoc", digits = 2)
outcome model AIC r2 (SD) RMSE (SD) term b [95% CI] SE t p
Body fat Null 712.18 Intercept 0.02 [-0.11, 0.14] 0.06 0.27 .791
Body fat ROI 712.92 0.02 (0.02) 1.00 (0.05) Intercept 0.02 [-0.10, 0.15] 0.06 0.38 .706
Body fat ROI 712.92 0.02 (0.02) 1.00 (0.05) Reward ROI -0.14 [-0.39, 0.11] 0.13 -1.12 .265
Body fat PEV 712.94 0.03 (0.04) 1.00 (0.07) Intercept 0.00 [-0.12, 0.13] 0.06 0.05 .960
Body fat PEV 712.94 0.03 (0.04) 1.00 (0.07) Craving PEV -0.09 [-0.26, 0.07] 0.08 -1.11 .268
Body fat Combined 714.49 0.03 (0.03) 1.01 (0.06) Intercept 0.01 [-0.12, 0.14] 0.07 0.19 .849
Body fat Combined 714.49 0.03 (0.03) 1.01 (0.06) Reward ROI -0.09 [-0.38, 0.19] 0.14 -0.66 .508
Body fat Combined 714.49 0.03 (0.03) 1.01 (0.06) Craving PEV -0.06 [-0.25, 0.13] 0.10 -0.65 .516

EMA enactment

enact_prop_fit = data.frame(model = c("ROI", "PEV", "COMBINED"),
                     r2 = c(mean(best_enact_prop_roi$resample$Rsquared), 
                            mean(best_enact_prop_pev$resample$Rsquared), 
                            mean(enact_prop_combined$resample$Rsquared)),
                     r2_sd = c(sd(best_enact_prop_roi$resample$Rsquared), 
                               sd(best_enact_prop_pev$resample$Rsquared), 
                               sd(enact_prop_combined$resample$Rsquared)),
                     RMSE = c(mean(best_enact_prop_roi$resample$RMSE), 
                              mean(best_enact_prop_pev$resample$RMSE), 
                              mean(enact_prop_combined$resample$RMSE)),
                     RMSE_sd = c(sd(best_enact_prop_roi$resample$RMSE), 
                                 sd(best_enact_prop_pev$resample$RMSE), 
                                 sd(enact_prop_combined$resample$RMSE)))

enact_table = broom::tidy(best_enact_prop_roi$finalModel, conf.int = TRUE) %>% 
  mutate(model = "ROI", 
         mod_sig = ifelse(pf(summary(best_enact_prop_roi$finalModel)$fstatistic[1], 
                             summary(best_enact_prop_roi$finalModel)$fstatistic[2], 
                             summary(best_enact_prop_roi$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", "")) %>%
  bind_rows(broom::tidy(best_enact_prop_pev$finalModel, conf.int = TRUE) %>% 
              mutate(model = "PEV",
                     mod_sig = ifelse(pf(summary(best_enact_prop_pev$finalModel)$fstatistic[1], 
                             summary(best_enact_prop_pev$finalModel)$fstatistic[2], 
                             summary(best_enact_prop_pev$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(enact_prop_combined$finalModel, conf.int = TRUE) %>% 
              mutate(model = "COMBINED",
                     mod_sig = ifelse(pf(summary(enact_prop_combined$finalModel)$fstatistic[1], 
                             summary(enact_prop_combined$finalModel)$fstatistic[2], 
                             summary(enact_prop_combined$finalModel)$fstatistic[3], lower = FALSE) < .05, "*", ""))) %>%
  bind_rows(broom::tidy(enact_prop_null, conf.int = TRUE) %>% mutate(model = "Null", mod_sig = "")) %>%
  rename("SE" = std.error,
         "t" = statistic,
         "p" = p.value) %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%.2f, %.2f]", estimate, conf.low, conf.high),
         p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p))),
         term = gsub("\\(Intercept\\)", "Intercept", term),
         term = gsub("_", " ", term),
         term = Hmisc::capitalize(term)) %>%
  left_join(., enact_aic) %>%
  left_join(., enact_prop_fit) %>%
  mutate(`r2 (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", r2, r2_sd), "--"),
         `RMSE (SD)` = ifelse(!model == "Null", sprintf("%.2f (%.2f)", RMSE, RMSE_sd), "--")) %>%
  unite(model, model, mod_sig, sep = "") %>%
  select(model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p) %>%
  arrange(AIC) %>%
  mutate(model = gsub("COMBINED", "Combined", model),
         outcome = "Enactment") %>%
    select(outcome, model, AIC, `r2 (SD)`, `RMSE (SD)`, term, `b [95% CI]`, SE, t, p)

enact_table %>%
  kable(format = "pandoc", digits = 2)
outcome model AIC r2 (SD) RMSE (SD) term b [95% CI] SE t p
Enactment PEV* 270.78 0.15 (0.13) 0.96 (0.16) Intercept -0.06 [-0.26, 0.13] 0.10 -0.64 .525
Enactment PEV* 270.78 0.15 (0.13) 0.96 (0.16) Reward PEV 0.37 [0.14, 0.59] 0.11 3.19 .002
Enactment Combined* 272.11 0.12 (0.11) 0.97 (0.11) Intercept -0.11 [-0.35, 0.12] 0.12 -0.97 .332
Enactment Combined* 272.11 0.12 (0.11) 0.97 (0.11) Value ROI 0.13 [-0.19, 0.45] 0.16 0.81 .419
Enactment Combined* 272.11 0.12 (0.11) 0.97 (0.11) Reward PEV 0.31 [0.05, 0.58] 0.13 2.33 .022
Enactment ROI* 275.56 0.11 (0.09) 0.99 (0.08) Intercept -0.16 [-0.39, 0.08] 0.12 -1.32 .189
Enactment ROI* 275.56 0.11 (0.09) 0.99 (0.08) Value ROI 0.32 [0.04, 0.60] 0.14 2.27 .026
Enactment Null 278.66 Intercept -0.01 [-0.21, 0.19] 0.10 -0.10 .918

combined

bind_rows(bmi_table, fat_table, enact_table) %>%
  kable(format = "pandoc", digits = 2)
outcome model AIC r2 (SD) RMSE (SD) term b [95% CI] SE t p
BMI Combined* 721.56 0.04 (0.05) 1.00 (0.13) Intercept -0.07 [-0.23, 0.08] 0.08 -0.94 .347
BMI Combined* 721.56 0.04 (0.05) 1.00 (0.13) Reward ROI -0.31 [-0.56, -0.05] 0.13 -2.37 .019
BMI Combined* 721.56 0.04 (0.05) 1.00 (0.13) Value ROI 0.12 [-0.08, 0.31] 0.10 1.17 .243
BMI Combined* 721.56 0.04 (0.05) 1.00 (0.13) Craving regulation PEV 0.14 [-0.05, 0.34] 0.10 1.47 .144
BMI ROI* 721.74 0.03 (0.05) 1.00 (0.11) Intercept -0.03 [-0.17, 0.11] 0.07 -0.43 .667
BMI ROI* 721.74 0.03 (0.05) 1.00 (0.11) Reward ROI -0.31 [-0.56, -0.05] 0.13 -2.35 .020
BMI ROI* 721.74 0.03 (0.05) 1.00 (0.11) Value ROI 0.16 [-0.03, 0.35] 0.10 1.61 .108
BMI PEV 723.45 0.02 (0.03) 0.99 (0.14) Intercept -0.05 [-0.19, 0.10] 0.07 -0.67 .507
BMI PEV 723.45 0.02 (0.03) 0.99 (0.14) Craving regulation PEV 0.16 [-0.03, 0.34] 0.09 1.64 .102
BMI Null 724.15 Intercept 0.01 [-0.11, 0.14] 0.06 0.21 .837
Body fat Null 712.18 Intercept 0.02 [-0.11, 0.14] 0.06 0.27 .791
Body fat ROI 712.92 0.02 (0.02) 1.00 (0.05) Intercept 0.02 [-0.10, 0.15] 0.06 0.38 .706
Body fat ROI 712.92 0.02 (0.02) 1.00 (0.05) Reward ROI -0.14 [-0.39, 0.11] 0.13 -1.12 .265
Body fat PEV 712.94 0.03 (0.04) 1.00 (0.07) Intercept 0.00 [-0.12, 0.13] 0.06 0.05 .960
Body fat PEV 712.94 0.03 (0.04) 1.00 (0.07) Craving PEV -0.09 [-0.26, 0.07] 0.08 -1.11 .268
Body fat Combined 714.49 0.03 (0.03) 1.01 (0.06) Intercept 0.01 [-0.12, 0.14] 0.07 0.19 .849
Body fat Combined 714.49 0.03 (0.03) 1.01 (0.06) Reward ROI -0.09 [-0.38, 0.19] 0.14 -0.66 .508
Body fat Combined 714.49 0.03 (0.03) 1.01 (0.06) Craving PEV -0.06 [-0.25, 0.13] 0.10 -0.65 .516
Enactment PEV* 270.78 0.15 (0.13) 0.96 (0.16) Intercept -0.06 [-0.26, 0.13] 0.10 -0.64 .525
Enactment PEV* 270.78 0.15 (0.13) 0.96 (0.16) Reward PEV 0.37 [0.14, 0.59] 0.11 3.19 .002
Enactment Combined* 272.11 0.12 (0.11) 0.97 (0.11) Intercept -0.11 [-0.35, 0.12] 0.12 -0.97 .332
Enactment Combined* 272.11 0.12 (0.11) 0.97 (0.11) Value ROI 0.13 [-0.19, 0.45] 0.16 0.81 .419
Enactment Combined* 272.11 0.12 (0.11) 0.97 (0.11) Reward PEV 0.31 [0.05, 0.58] 0.13 2.33 .022
Enactment ROI* 275.56 0.11 (0.09) 0.99 (0.08) Intercept -0.16 [-0.39, 0.08] 0.12 -1.32 .189
Enactment ROI* 275.56 0.11 (0.09) 0.99 (0.08) Value ROI 0.32 [0.04, 0.60] 0.14 2.27 .026
Enactment Null 278.66 Intercept -0.01 [-0.21, 0.19] 0.10 -0.10 .918

VIF

# BMI
car::vif(bmi_combined$finalModel)
##             reward_ROI              value_ROI craving_regulation_PEV 
##               1.101738               1.185011               1.085275
# %fat
car::vif(fat_combined$finalModel)
##  reward_ROI craving_PEV 
##    1.303633    1.303633
# enactment
car::vif(enact_prop_combined$finalModel)
##  value_ROI reward_PEV 
##   1.348297   1.348297